
# Required packages
library(mfx)
library(ggplot2)
library(gridExtra)
library(reshape)
library(maps)
library(interactions)
library(brglm)
library(miceadds)
library(dplyr)
library(stats)
library(Zelig)

setwd()

#######################################################################################
## Function to create tables
#######################################################################################

# create preamble with number of columns corresponding to number of models
table_me <- function(formatted.results, num.obs, caption, digits){
  m <- length(formatted.results)
  c <- paste(rep("D{.}{.}{2}", m), collapse = '')
  
  preamble <- paste('\\begin{table}[!ht]
                    \\caption{', caption, '}
                    \\label{} 
                    \\begin{center}
                    \\begin{tabular}{ l ', c, ' } 
                    \\hline \\hline
                    ', sep = '')
  
  # row with column headers
  multicolumns <- "& \\multicolumn{ 1 }{ c }{ Model 1 }"
  for(i in 2:m){
    multicolumns <-  paste(multicolumns, "& \\multicolumn{ 1 }{ c }{ Model", i, "} ", collapse = '')
  }
  multicolumns <- paste(multicolumns, '\\\\\\hline', collapse = '')
  
  # paste together
  header <- paste(preamble, multicolumns, collapse = '')
  
  # create list of unique variables
  vars <- row.names(formatted.results[[1]])
  for(j in 2:m){
    vars <- unique(c(vars, row.names(formatted.results[[j]])))
  }
  
  # loop to create lines for each var
  table <- header
  for(k in 1:length(vars)){
    line1 <- paste(vars[k], '&', collapse = '')
    line2 <- paste('    ', '&', collapse = '')
    # loop to extract coefs for 1 through m-1
    for(j in 1:(m-1)){
      # if to determine if coef is present in model m
      if(vars[k] %in% row.names(formatted.results[[j]])){
        # if so paste coef and se
        line1 <- paste(line1, format(round(formatted.results[[j]][row.names(formatted.results[[j]]) == vars[k], 1], digits), nsmall = digits))
        line2 <- paste(line2, '(', format(round(formatted.results[[j]][row.names(formatted.results[[j]]) == vars[k], 2], digits), nsmall = digits), ') &')
        # if to determine significance
        if(formatted.results[[j]][row.names(formatted.results[[j]]) == vars[k], 3] < 0.05){
          # if so add star
          line1 <- paste(line1, '^* &', sep = '')
        }else{
          line1 <- paste(line1, '&')
        }
      }else{
        line1 <- paste(line1, ' & ')
        line2 <- paste(line2, '    &')
        
      }
    } 
    # add coefs for last model
    # if to determine if coef is present in model m
    if(vars[k] %in% row.names(formatted.results[[m]])){
      # if so paste coef and se
      line1 <- paste(line1, format(round(formatted.results[[m]][row.names(formatted.results[[m]]) == vars[k], 1], digits), nsmall = digits))
      line2 <- paste(line2, '(', format(round(formatted.results[[m]][row.names(formatted.results[[m]]) == vars[k], 2], digits), nsmall = digits), ')')
      # if to determine significance
      if(formatted.results[[m]][row.names(formatted.results[[m]]) == vars[k], 3] < 0.05){
        # if so add star
        line1 <- paste(line1, '^* ', sep = '')
      }
    }
    
    # add final line breaks
    line1 <- paste(line1, '\\\\')
    line2 <- paste(line2, '\\\\')
    table <- paste(table, line1, line2, sep = '\n')
  }
  # add number of observations
  obs <- paste('$N$', num.obs[[1]], sep = ' & ')
  for(l in 2:m){
    obs <- paste(obs, num.obs[[l]], sep = ' & ')
  }
  obs <- paste(obs, '\\\\')
  table <- paste(table, '\\hline', '\n', obs, '\\hline \\hline')
  table <- paste(table, '\\end{tabular}', '\\end{center}', '\\caption*{\\scriptsize{}}',
                 '\\end{table}', sep = '\n')
  return(writeLines(table))
}


#######################################################################################
## Main Text
#######################################################################################

#######################################################################################
## Figure 1 
#######################################################################################

# Load data
bon <- read.csv('bod_fortune_500_DIME_cont_records.csv')
bon <- bon[bon$cycle == '2010', ]

# candidates on house energy and commerce

energy_commerce <- c('BALDWIN, TAMMY', 'BARROW, JOHN J', 'BARTON, JOE LINUS', 'BLACKBURN, MARSHA W',
                     'BLUNT, ROY', 'BONO, MARY', 'BOUCHER, FREDERICK C "RICK"', 'BRALEY, BRUCE L.',
                     'BURGESS, MICHAEL CLIFTON', 'BUTTERFIELD, G. K.', 'BUYER, STEPHEN E', 'CAPPS, LOIS',
                     'CASTOR, KATHY', 'CHRISTENSEN, DONNA M', 'DEGETTE, DIANA LOUISE', 'DINGELL, JOHN D.',
                     'DOYLE, MICHAEL F JR', 'ENGEL, ELIOT LANCE', 'ESHOO, ANNA', 'GINGREY, J PHILLIP MD',
                     'GONZALEZ, CHARLES A', 'GORDON, BARTON J', 'GREEN, GENE', 'GRIFFITH, PARKER REP.',
                     'HALL, RALPH MOODY', 'HARMAN, JANE', 'HILL, BARON PAUL', 'INSLEE, JAY R', 'LATTA, ROBERT',
                     'MARKEY, EDWARD J', 'MATHESON, JAMES DAVID', 'MATSUI, DORIS', 'MCNERNEY, JERRY',
                     'MELANCON, CHARLES J', 'MURPHY, TIM', 'MURPHY, CHRISTOPHER SCOTT', 'MYRICK, SUE',
                     'PALLONE, FRANK JR', 'PITTS, JOSEPH RUSSELL', 'RADANOVICH, GEORGE', 'ROGERS, MIKE',
                     'ROSS, MICHAEL AVERY', 'RUSH, BOBBY L', 'SARBANES, JOHN PETER SPYROS', 'SCALISE, STEPHEN J', 'SCHAKOWSKY, JANICE D',
                     'SHADEGG, JOHN BARDEN', 'SHIMKUS, JOHN MONDY', 'SPACE, ZACHARY T', 'STEARNS, CLIFFORD B', 'STUPAK, BART',
                     'SULLIVAN, JOHN', 'SUTTON, BETTY S MS.', 'TERRY, LEE R', 'UPTON, FREDERICK STEPHEN', 'WAXMAN, HENRY A.',
                     'WEINER, ANTHONY DAVID', 'WELCH, PETER', 'WHITFIELD, EDWARD')

bon <- bon[bon$cycle == 2010 & bon$recipient.name %in% energy_commerce, ]

ids <- read.csv('bod_fortune_500_DIME.csv')[, c("dime.cid", "corp.name", "sector", "industry")]
names(ids) <- c("dime.cid", "corpname", "sector", "industry")

bon <- merge(bon, ids, by = c("dime.cid", "corpname"), all.x = T)
giving <- bon %>% group_by(corpname, recipient.party) %>% summarize(Giving = sum(amount, na.rm = T))
totalgiving <- bon %>% group_by(corpname) %>% summarize(TotalGiving = sum(amount, na.rm = T))

totalgiving <- merge(totalgiving, giving[giving$recipient.party == 200, c("corpname", "Giving")], all.x = T)
totalgiving$Giving <- ifelse(is.na(totalgiving$Giving), 0, totalgiving$Giving)
totalgiving$PropDem <- totalgiving$Giving / totalgiving$TotalGiving

totalgiving <- merge(totalgiving, na.omit(unique(bon[, c("corpname", "sector", "industry")])))
names(totalgiving)[2] <- 'TotalContributions'

means <- totalgiving %>% group_by(sector) %>% summarize(dem = sum(Giving, na.rm = T), total = sum(TotalContributions, na.rm = T))
means$Mean <- means$dem / means$total
totalgiving <- merge(totalgiving, means[, c("sector", "Mean")], by = "sector", all.x = T)
totalgiving$sector <- factor(totalgiving$sector, levels = unique(totalgiving$sector[order(totalgiving$Mean)]))

set.seed(5)
ggplot(totalgiving[totalgiving[,3] > 0, ], aes(sector, PropDem, size = TotalContributions, color = PropDem)) + geom_point() + geom_jitter(width = 0.00, height = 0.05) +
  scale_size(range = c(1,6)) + scale_colour_gradient(low = "#D73027", high = "#4575B4", guide = 'none') + theme(text = element_text(size=20), legend.title=element_text(size=16)) +
  geom_hline(yintercept = mean(giving$PropDem, na.rm = T), col = "black", lty = 3) + coord_flip() +
  labs(x = "Sector", y = "Proportion Contributions to Democrats", size = 'Total Contributions') 




#######################################################################################
## Figure 2 
#######################################################################################

# TOP PANEL
osiris <- read.csv("osirisdata-3-29-17.CSV", header = T)
dups <- osiris[duplicated(osiris$OS_ID_NUMBER) == T, ]
osirisnodups <- osiris[!duplicated(osiris$OS_ID_NUMBER), ] 
osirisnodups$PlantMach <- as.numeric(as.character(osirisnodups$PlantMach))
osirisnodups$NetPPE <- as.numeric(as.character(osirisnodups$NetPPE))
ppe <- osirisnodups %>% group_by(CAICS12COD) %>% summarize(MeanNetPPE = mean(NetPPE, na.rm = T),
                                                           SDNetPPE = sd(NetPPE, na.rm = T),
                                                           MeanPlantMach = mean(PlantMach, na.rm = T),
                                                           SDPlantMach = sd(PlantMach, na.rm = T),
                                                           NoFirms = n())

ppe$Manufacturing <- ifelse(substr(ppe$CAICS12COD, 1, 2) == '31' | substr(ppe$CAICS12COD, 1, 2) == '32' | substr(ppe$CAICS12COD, 1, 2) == '33', 'Manufacturing', 'Non-Manufacturing')
ppe$Manufacturing <- factor(ppe$Manufacturing, levels = c('Manufacturing', 'Non-Manufacturing'))
ppe <- ppe[!is.na(ppe$Manufacturing), ]


set.seed(5)
plot1 <- ggplot(data=ppe, aes(log(MeanNetPPE), log(SDNetPPE), size = NoFirms, col = Manufacturing)) + 
  geom_point()  + geom_jitter(width = 2, height = 2) +
  scale_x_continuous(labels = scales::comma) + scale_y_continuous(labels = scales::comma) +
  labs(x="(Log) Mean Net PP&E by Industry", y="(Log) St. Dev. Net PP&E by Ind.") +
  scale_colour_manual(values = c('red', "gray")) +
  theme(legend.position = c(.15, 0.75), text = element_text(size=15)) + guides(color = FALSE, size=guide_legend("Industry Size"))


set.seed(5)
plot2 <- ggplot(data=ppe, aes(log(MeanPlantMach), log(SDPlantMach), size = NoFirms, col = Manufacturing)) + 
  geom_point()  + geom_jitter(width = 2, height = 2) +
  scale_x_continuous(labels = scales::comma) + scale_y_continuous(labels = scales::comma) +
  labs(x="(Log) Mean Plant & Machinery by Industry", y="(Log) St. Dev. Plant & Mach. by Ind.") +
  scale_colour_manual(values = c('red', "gray")) +
  theme(legend.position = c(.8, 0.2), text = element_text(size=15)) + guides(size = FALSE, color=guide_legend("Industry Type"))

# BOTTOM PANEL
mecs1 <- read.csv('Table 8.2-Table 1.csv')[8:98,]

technologies <- t(as.matrix(mecs1))[,1]
technologies <- unique(technologies)[2:6]
technologies[1] <- substr(technologies[1], 1, 44)
technologies[2] <- substr(technologies[2], 1, 61)


mecs1.a <- mecs1[16:91, 1:6]
mecs1.b <- mecs1[16:91, c(1:3, 7:9)]
mecs1.c <- mecs1[16:91, c(1:3, 10:12)]
mecs1.d <- mecs1[16:91, c(1:3, 13:15)]
mecs1.e <- mecs1[16:91, c(1:3, 16:18)]

names(mecs1.a) <- c('NAICS', 'Industry', 'Establishments', 'TechAInUse', 'TechANotInUse', "TechADK")
names(mecs1.b) <- c('NAICS', 'Industry', 'Establishments', 'TechBInUse', 'TechBNotInUse', "TechBDK")
names(mecs1.c) <- c('NAICS', 'Industry', 'Establishments', 'TechCInUse', 'TechCNotInUse', "TechCDK")
names(mecs1.d) <- c('NAICS', 'Industry', 'Establishments', 'TechDInUse', 'TechDNotInUse', "TechDDK")
names(mecs1.e) <- c('NAICS', 'Industry', 'Establishments', 'TechEInUse', 'TechENotInUse', "TechEDK")

mecs1 <- merge(mecs1.a, mecs1.b, by = c('NAICS', 'Industry', 'Establishments'))
mecs1 <- merge(mecs1, mecs1.c, by = c('NAICS', 'Industry', 'Establishments'))
mecs1 <- merge(mecs1, mecs1.d, by = c('NAICS', 'Industry', 'Establishments'))
mecs1 <- merge(mecs1, mecs1.e, by = c('NAICS', 'Industry', 'Establishments'))

rm(mecs1.a, mecs1.b, mecs1.c, mecs1.d, mecs1.e)

for(i in 3:18){
  mecs1[, i] <- as.numeric(gsub(',', '', as.character(mecs1[, i])))
}

mecs1$PropA <- mecs1$TechAInUse / (mecs1$TechAInUse + mecs1$TechANotInUse)
mecs1$PropB <- mecs1$TechBInUse / (mecs1$TechBInUse + mecs1$TechBNotInUse)
mecs1$PropC <- mecs1$TechCInUse / (mecs1$TechCInUse + mecs1$TechCNotInUse)
mecs1$PropD <- mecs1$TechDInUse / (mecs1$TechDInUse + mecs1$TechDNotInUse)
mecs1$PropE <- mecs1$TechEInUse / (mecs1$TechEInUse + mecs1$TechENotInUse)

mecs1 <- mecs1[nchar(trimws(as.character(mecs1$NAICS))) == 6, ]

mecs1$`Computer Control of Building Wide Environment` <- mecs1$PropA
mecs1$`Computer Control of Processes or Major Energy-Using Equipment` <- mecs1$PropB
mecs1$`Waste Heat Recovery` <- mecs1$PropC
mecs1$`Adjustable Speed Motors` <- mecs1$PropD
mecs1$`Oxy Fuel Firing` <- mecs1$PropE

mecs1.melt <- melt(mecs1[, c(1, 19:23)])
names(mecs1.melt) <- c("NAICS", "variable", "Proportion")

technames <- as.data.frame(cbind(c('PropA', 'PropB', 'PropC', 'PropD', 'PropE'), 
                                 c('Automated Building Environment', 'Automated Processes and Equipment',
                                   'Waste Heat Recovery', 'Adjustment Speed Motors', 'Oxy Fuel Firing')))
names(technames) <- c("variable", "Technology Employed")
mecs1.melt <- merge(mecs1.melt, technames, by = 'variable')


plot3 <- ggplot(mecs1.melt, aes(Proportion, group = `Technology Employed`, fill = `Technology Employed`)) + 
  geom_density(position="stack") +
  theme(legend.title = element_blank(), text = element_text(size=15), legend.position = c(0.7, 0.75)) + 
  labs(x = 'Use of Energy Saving Technologies', y = 'Density') +
  scale_fill_grey() + xlim(c(0,1)) 

mecs2 <- read.csv('Table 10.7-Table 1.csv')[17:136, 1:5]
mecs3 <- read.csv('Table 10.11-Table 1.csv')[17:136, 1:5]

names(mecs2) <- c('NAICS', 'Industry', 'Establishments', 'Switchable', 'NotSwitchable')
names(mecs3) <- c('NAICS', 'Industry', 'Establishments', 'Switchable', 'NotSwitchable')

mecs2 <- mecs2[nchar(trimws(as.character(mecs2$NAICS))) == 6, ]
mecs3 <- mecs3[nchar(trimws(as.character(mecs3$NAICS))) == 6, ]

for(i in 4:5){
  mecs2[, i] <- as.numeric(gsub(',', '', as.character(mecs2[, i]))) 
  mecs3[, i] <- as.numeric(gsub(',', '', as.character(mecs3[, i]))) 
}

mecs2$Electricity_Switch <- mecs2$Switchable / (mecs2$Switchable + mecs2$NotSwitchable)
mecs3$Coal_Switch <- mecs3$Switchable / (mecs3$Switchable + mecs3$NotSwitchable)

mecs4 <- merge(mecs2[, c(1,2, 6)], mecs3[, c(1,2,3, 6)], by = c('NAICS', 'Industry'))
mecs4$Establishments <- as.numeric(as.character(mecs4$Establishments))

mecs4$NAICS2 <- substr(trimws(as.character(mecs4$NAICS)), 1, 2)

set.seed(5)
plot4 <- ggplot(mecs4, aes(Electricity_Switch, Coal_Switch, size = Establishments)) + geom_point() +
  geom_jitter(width = 0.075, height = 0.075) + 
  labs(x = 'Electricity to Alternative Energy', y = 'Coal to Alternative Energy', size = 'Industry Size') +
  theme(text = element_text(size=15), legend.position = c(0.85, 0.25))

grid.arrange(plot1, plot2, plot3, plot4, ncol = 2)

#######################################################################################
## Figure 3 
#######################################################################################

# LEFT PANEL
fulldata <- read.csv('Enemies_Replication.csv')

statedata <- na.omit(unique(fulldata[, c("State", "StateName", "StateEnergy")]))

names(statedata) <- c("state", "region", "Emissions")
statedata$region <- as.character(statedata$region)
statedata$region <- tolower(statedata$region)

states <- map_data("state")
tfmerged <- merge(statedata, states, sort = FALSE, by = "region")
tfmerged <- tfmerged[order(tfmerged$order), ]

plot <- qplot(long, lat, data = tfmerged, group = group, fill = -log(Emissions), 
              geom="polygon", ylab = "", xlab = "")
plot1 <- plot + theme_bw() + theme(
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks = element_blank()) +
  labs(x="State-Level Carbon Emissions per Capita", y="") + 
  theme(legend.position="none", axis.title.x = element_text(margin = margin(t = 16, r = 0, b = 0, l = 0)), text = element_text(size=20))

# RIGHT PANEL

plot2 <- ggplot(data=fulldata, aes(StateEnergy)) + 
  geom_density(alpha=.8, fill = "grey") + labs(x="", y="Firm Weighted Density") + geom_vline(xintercept = 10.62477, linetype = "longdash", color = "black") + 
  geom_vline(xintercept = 24.46256, linetype = "longdash", color = "black") + scale_x_continuous('Carbon Emissions per Capita', limits = c(0, 50)) + 
  geom_text(data = NULL, x = 17, y = 0.125, label = "Connecticut", size=8) + 
  scale_size(1) + geom_text(data = NULL,  x = 28, y = 0.10, label = "Texas", size=8) + theme_bw() + theme(text = element_text(size=20))

grid.arrange(plot1, plot2, ncol = 2)

#######################################################################################
## Figure 4 
#######################################################################################

# load stock returns and subset to active
stocks <- read.csv('585fed133b24cb5c.csv')
stocks <- stocks[stocks$SECSTAT == 'R', ]   

# reformat variables, create index and time vars
stocks$COMNAM <- as.character(stocks$COMNAM)
stocks$COMNAM <- ifelse(stocks$COMNAM == 'DOW 30 PREMIUM & DIV INC FD INC', 'Dow', stocks$COMNAM)
stocks$date <- as.Date(stocks$date, format = '%m/%d/%Y')


# calculate returns
stocks$RETX <- (stocks$PRC - stocks$OPENPRC) / stocks$OPENPRC

dow <- stocks[stocks$COMNAM == 'Dow', c('RETX', "date")]
names(dow) <- c('Index', 'date')
stocks <- merge(stocks, dow, by = 'date', all.x = T)
stocks$t <- stocks$date - as.Date('2009-06-29')

# set up event and estimation windows
est.win <- c(-30, -6)
ev.win <- c(0, 4)

# define function to calculate cumulative abnormal returns
car.function <- function(est.win, ev.win, firm, df){
  # subset to country i in the estimation window
  est.sub <- df[df$t >= est.win[1] & df$t <= est.win[2] & df$COMNAM == firm, ]
  if(nrow(est.sub) > 5){
    
    # regress spread on index + covariates if applicable
      fit <- lm(RETX ~ Index, data = est.sub) 
    # if any coefficients not well defined then say so
    if(sum(is.na(fit$coefficients)) > 0){print(paste("Cannot estimate ar for", firm, sep = " "))}
    success <- ifelse(sum(!is.na(fit$coefficients)) > 0, 1, 0)
    
    # condition on model fit to eliminate those that aren't well-predicted by index
    if(success == 1){
        # subset to event window + estimation window for country i
        ev.sub <- df[((df$t >= est.win[1] & df$t <= est.win[2]) | (df$t >= ev.win[1] & df$t <= ev.win[2])) & 
                       df$COMNAM == firm, ]
        ev.sub$Pred <- as.numeric(as.matrix(cbind(1, ev.sub$Index)) %*% fit$coefficients)
          }  
        
        # calculate abnormal return
        ev.sub$ar <- ev.sub$RETX - ev.sub$Pred
        
        # calculate standard deviation of abnormal return for estimation window only
        s.ar <- sqrt((1/(est.win[2] - est.win[1]))*sum((ev.sub$ar[ev.sub$t >= est.win[1] & ev.sub$t <= est.win[2]])^2, na.rm = T))
        
        # cumulative abnormal return
        car <- sum(ev.sub$ar[ev.sub$t > est.win[2]])
        
        # standard deviation of cumulative abnormal return for estimation window
        s.car <- sqrt(ev.win[2] - ev.win[1] + 1)*s.ar
        
        # calculate significance of cumulative abnormal returns for country i
        car.sub <- c(as.character(firm), car, s.car, round(car/s.car, digits = 2), 
                     round(pt(-1*abs(car/s.car), ev.win[2] - ev.win[1]), digits = 2), 
                     ifelse(pt(-1*abs(car/s.car), ev.win[2] - ev.win[1]) <= 0.05, '*', ''),
                     summary(fit)$adj.r.squared)
        
        return(car.sub)
  } else {print("no obs in window")}}

# loop to calculate car for all firms
car.holder <- rep(NA, 7)
for(i in 1:length(unique(stocks$COMNAM))){
  unique(stocks$COMNAM)[i]
  car <- car.function(est.win, ev.win, unique(stocks$COMNAM)[i], stocks)
  car.holder <-   rbind(car.holder, car)
} 


# joint significance of abnormal returns
caar.sub <- NULL
# test significance of caar
car.holder <- as.data.frame(car.holder)
car.holder <- car.holder[!is.na(car.holder$V2), ]
car.holder[,2] <- as.numeric(as.character(car.holder[,2]))
car.holder[,3] <- as.numeric(as.character(car.holder[,3]))

caar <- mean(car.holder3[,2], na.rm = T)
s.caar <- (1/nrow(car.holder3)^2)*sqrt(sum((car.holder3[,3])^2, na.rm = T)) # see http://web.worldbank.org/archive/website01004/WEB/IMAGES/12852832.PDF
caar.sub <- c(caar, s.caar, caar/s.caar, 
              pt(-abs(caar/s.caar), ev.win[2] - ev.win[1]),
              ifelse(pt(-abs(caar/s.caar), ev.win[2] - ev.win[1]) <= 0.05, '*', ''))
names(caar.sub) <- c("caar", "st. error", "z-stat", "p-value", "significance")


# merge abnormal returns with industry data

industries <- unique(stocks[, c("PERMNO", "COMNAM", "TSYMBOL", "NAICS")])
cars <- car.holder[, c(1,2, 6)]
names(cars) <- c("COMNAM", "car3", "sig3")
cars$Sig <- ifelse(cars$sig3 == '*', 'p < 0.05', 'p > 0.05')
cars$Sig <- factor(cars$Sig, levels = c('p < 0.05', 'p > 0.05'))
df <- merge(industries, cars, by = "COMNAM", all.x = T)    


# add naics codes
codes <- read.csv('6-digit_2012_Codes.csv')[, 1:2]
names(codes) <- c("NAICS", "NAICS_NAME")
df <- merge(df, codes, by = 'NAICS', all.x = T)                                            

# 2-6 digit codes
codes26 <- read.csv('2-6 digit_2017_Codes.csv')[,2:3]
names(codes26) <- c("NAICS", "NAICS_NAME")

codes2 <- codes26[nchar(as.character(codes26$NAICS)) == 2, ]
codes4 <- codes26[nchar(as.character(codes26$NAICS)) == 4, ]

df$NAICS2 <- substr(as.character(df$NAICS), 1, 2)
df$NAICS4 <- substr(as.character(df$NAICS), 1, 4)

names(codes2) <- c("NAICS2", "NAICS2_NAME")
names(codes4) <- c("NAICS4", "NAICS4_NAME")

df <- merge(df, codes2, by = "NAICS2", all.x = T)
df <- merge(df, codes4, by = "NAICS4", all.x = T)


# random reformatting
df <- df[!is.na(df$car3) & !is.na(df$sig3), ]
df <- df[!is.na(df$NAICS), ]
df$NAICS4_NAMEB <- gsub('Manufacturing', '', df$NAICS4_NAME)

# summarize by industry group
sum <- df %>% 
  group_by(NAICS4_NAMEB) %>% 
  summarize(NAICS4 = max(NAICS4), mean = mean(car3, na.rm = T), sd = sd(car3, na.rm = T), NAICS2 = max(NAICS2))

df_manuf <- df[df$NAICS2 %in% 31:33, ]
naics4_list <- unique(df_manuf$NAICS4)

# create plots
plot1 <- ggplot(df_manuf[df_manuf$NAICS4 %in% naics4_list[1:28], ], aes(x = NAICS4, y = car3, col = sig3), size = .75) + 
  geom_point(size=2) + geom_hline(yintercept = 0, lty = 2, col = 'black') +
  scale_colour_manual(values = c("grey", "black")) +
  theme(legend.position="none", text = element_text(size=20)) +
  labs(x = "NAICS4 Industry", y = 'Abnormal Return (%)') + coord_flip() + ylim(-0.5, 0.5) +
  geom_point(aes(x = NAICS4, y = mean), data = sum[sum$NAICS4 %in% naics4_list[1:28], ], color="red", size = 2, shape = 17)

plot2 <- ggplot(df_manuf[df_manuf$NAICS4 %in% naics4_list[29:56], ], aes(x = NAICS4, y = car3, col = sig3), size = .75) + 
  geom_point(size=2) + geom_hline(yintercept = 0, lty = 2, col = 'black') +
  scale_colour_manual(values = c("grey", "black")) +
  theme(legend.position="none", text = element_text(size=20)) +
  labs(x = "NAICS4 Industry", y = 'Abnormal Return (%)') + coord_flip() + ylim(-0.5, 0.5) +
  geom_point(aes(x = NAICS4, y = mean), data = sum[sum$NAICS4 %in% naics4_list[29:56], ], color="red", size = 2, shape = 17)

plot3 <- ggplot(df_manuf[df_manuf$NAICS4 %in% naics4_list[57:83], ], aes(x = NAICS4, y = car3, col = Sig), size = .75) + 
  geom_point(size=2) + geom_hline(yintercept = 0, lty = 2, col = 'black') +
  scale_colour_manual(values = c("black", "grey")) +
  theme(legend.title = element_blank(), legend.position = c(0.8, 0.08), text = element_text(size=20)) +
  labs(x = "NAICS4 Industry", y = 'Abnormal Return (%)') + coord_flip() + ylim(-0.5, 0.5) +
  geom_point(aes(x = NAICS4, y = mean), data = sum[sum$NAICS4 %in% naics4_list[57:83], ], color="red", size = 2, shape = 17)

grid.arrange(plot1, plot2, plot3, ncol = 3)


#######################################################################################
## Figure 6 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# TOP PANEL

# Model
fit1 <- logitmfx(Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                   NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, data = fulldata, atmean = T, clustervar1 = "CAICS12COD")
mfx <- as.data.frame(fit1$mfxest)

z_score <- qnorm(0.975)
mfx$UpperCI <- mfx[,1] + mfx[,2] * z_score
mfx$LowerCI <- mfx[,1] - mfx[,2] * z_score
mfx$Variable <- row.names(mfx)
names(mfx) <- c("MarginalEffect", "StError", "z", "P", "UpperCI", "LowerCI", "Variable")

# Re-name variables
mfx[1,7] <- "Competitor Costs (Coal Reserves)"
mfx[2,7] <- "Own Costs (Coal Reserves)"
mfx[5,7] <- "PPE"
mfx[7,7] <- "Energy Intensive Manufacturing"
mfx[8,7] <- "Multinational"
mfx[9,7] <- "Electricity Provider"

# Re-order factors
mfx[,7] <- factor(mfx[,7], levels = c(
  "PPE",
  "MarketCap",
  "Productivity",
  "MarketShare",
  "Electricity Provider",
  "Energy Intensive Manufacturing",
  "Multinational",
  "Own Costs (Coal Reserves)",
  "Competitor Costs (Coal Reserves)"
))

# Set colors
myred <- "#D73027"


# coefplot
plot1 <- ggplot(mfx[c(1,2,8,3),], aes(Variable, MarginalEffect)) +
  geom_errorbar(aes(x = Variable, ymin = LowerCI, ymax = UpperCI), size = 2, width = .5, color = myred) + 
  scale_x_discrete('') +
  scale_y_continuous('') +
  geom_point(aes(x = Variable, y = MarginalEffect), size = 5, color = myred) +
  geom_hline(yintercept = 0) + 
  theme_bw() + theme(text = element_text(size=20)) + 
  coord_flip() 

# BOTTOM PANEL

# Model
fit1 <- logitmfx(Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + NetPPE + 
                   MarketCap + HighManufEnergy + FSubs + electricpowergeneration, 
                 data = fulldata, atmean = T, clustervar1 = "CAICS12COD")
mfx <- as.data.frame(fit1$mfxest)

z_score <- qnorm(0.975)
mfx$UpperCI <- mfx[,1] + mfx[,2] * z_score
mfx$LowerCI <- mfx[,1] - mfx[,2] * z_score
mfx$Variable <- row.names(mfx)
names(mfx) <- c("MarginalEffect", "StError", "z", "P", "UpperCI", "LowerCI", "Variable")

# Re-name variables
mfx[1,7] <- "Competitor Costs (Coal Reliance)"
mfx[2,7] <- "Own Costs (Coal Reliance)"
mfx[5,7] <- "PPE"
mfx[7,7] <- "Energy Intensive Manufacturing"
mfx[8,7] <- "Multinational"
mfx[9,7] <- "ElectricityProvider"

# Re-order factors
mfx[,7] <- factor(mfx[,7], levels = c(
  "PPE",
  "MarketCap",
  "Productivity",
  "MarketShare",
  "ElectricityProvider",
  "Energy Intensive Manufacturing",
  "Multinational",
  "Own Costs (Coal Reliance)",
  "Competitor Costs (Coal Reliance)"
))


# Set colors
myblue <- "#4575B4"


# coefplot
plot2 <- ggplot(mfx[c(1,2,8,3),], aes(Variable, MarginalEffect)) +
  scale_x_discrete('') +
  scale_y_continuous('', limits = c(-.15, .35)) +
  theme_bw() + 
  geom_errorbar(aes(x = Variable, ymin = LowerCI, ymax = UpperCI), size = 2, width = .5, color = myblue) + 
  geom_point(aes(x = Variable, y = MarginalEffect), size = 5, color = myblue) +
  geom_hline(yintercept = 0) + 
  theme(text = element_text(size=20)) +
  coord_flip() 

grid.arrange(plot1, plot2, ncol = 1)



#######################################################################################
## Figure 7
#######################################################################################

# read data
p_small <- read.csv('p_small.csv')
s_small <- read.csv('s_small.csv')

# load maps
countries <- map_data("world")
tfmerged <- merge(p_small, countries, sort = FALSE, by = "region", all.y = T)
tfmerged <- tfmerged[order(tfmerged$order), ]

plot <- qplot(long, lat, data = tfmerged, group = group, fill = ClimFinP, 
              geom="polygon", ylab = "", xlab = "")
plot1 <- plot + theme_bw() + theme(
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks = element_blank()) +
  scale_fill_gradient(name = "Millions USD") +
  labs(x="Climate Finance (Principal Objective)", y="", legend = 'Millions USD') + 
  theme(axis.title.x = element_text(margin = margin(t = 16, r = 0, b = 0, l = 0)),
        text = element_text(size=20))

# RIGHT PANEL
tfmerged <- merge(s_small, countries, sort = FALSE, by = "region", all.y = T)
tfmerged <- tfmerged[order(tfmerged$order), ]

plot <- qplot(long, lat, data = tfmerged, group = group, fill = ClimFinS, 
                geom="polygon", ylab = "", xlab = "")
plot2 <- plot + theme_bw() + theme(
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks = element_blank()) +
  scale_fill_gradient(name = "Millions USD") +
  labs(x="Climate Finance (Significant Objective)", y="", legend = 'Millions USD') + 
  theme(axis.title.x = element_text(margin = margin(t = 16, r = 0, b = 0, l = 0)),
        text = element_text(size=20))

grid.arrange(plot1, plot2, ncol = 2)

#######################################################################################
## Figure 8 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# TOP PANEL

probs <- quantile(fulldata$LogPenetration, probs = c(.5, .6, .7, .8, .9), na.rm = T)

plot1 <- interact_plot(glm(formula = Support ~ CompShareZipCoal*LogPenetration + ZipCoalField + MarketShare + Productivity + 
                             NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, data = fulldata, family = "binomial"),
                       pred = "CompShareZipCoal", modx = "LogPenetration", modx.values = probs, 
                       modx.labels = c('50th', '60th', '70th', '80th', '90th'),
                       legend.main = 'Log Penetration \n(Percentile)', x.label = 'Competitor Costs (Coal Reserves)')

plot2 <- interact_plot(glm(formula = Support ~ CompShareZipCoal*teii + ZipCoalField + MarketShare + Productivity + 
                             NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, data = fulldata, family = "binomial"),
                       pred = "CompShareZipCoal", modx = "teii", legend.main = 'Trade Exposed, \nEnergy Intensive',
                       x.label = 'Competitor Costs (Coal Reserves)', colors = 'Blues')



# exports

probs <- quantile(fulldata$LogExports, probs = c(.5, .6, .7, .8, .9), na.rm = T)

plot3 <- interact_plot(glm(formula = Support ~ CompShareZipCoal*LogExports + ZipCoalField + MarketShare + Productivity + 
                             NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, data = fulldata, family = "binomial"),
                       pred = "CompShareZipCoal", modx = "LogExports", modx.values = probs,
                       modx.labels = c('50th', '60th', '70th', '80th', '90th'),
                       legend.main = 'Log Exports \n(Percentile)', x.label = 'Competitor Costs (Coal Reserves)', colors = 'Reds')

plot4 <- interact_plot(glm(formula = Support ~ CompShareZipCarbon*LogExports + ZipCoalField + MarketShare + Productivity + 
                             NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, data = fulldata, family = "binomial"),
                       pred = "CompShareZipCarbon", modx = "LogExports",  modx.values = probs, 
                       modx.labels = c('50th', '60th', '70th', '80th', '90th'),
                       legend.main = 'Log Exports \n(Percentile)',
                       x.label = 'Competitor Costs (Carbon Reliance)', colors = 'Reds')

grid.arrange(plot1, plot2, plot3, plot4, ncol = 2)



#######################################################################################
## Figure 9 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

vars <- c('FSubs', 'FinalGood', 'CompShareZipCoal', 'ZipCoalField', 'MarketShare', 'Productivity', 'NetPPE', 'MarketCap',
          'HighManufEnergy', 'electricpowergeneration', 'Support')
dat <- fulldata[complete.cases(fulldata[, vars]), ]

dat$MNC1Final0 <- ifelse(dat$FSubs == 1 & dat$FinalGood == 0, 1, 0)
dat$MNC0Final1 <- ifelse(dat$FSubs == 0 & dat$FinalGood == 1, 1, 0)
dat$MNC1Final1 <- ifelse(dat$FSubs == 1 & dat$FinalGood == 1, 1, 0)
formula <- 'Support ~ MNC1Final0 + MNC0Final1 + MNC1Final1 + CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
NetPPE + MarketCap + HighManufEnergy + electricpowergeneration'

fit1 <- brglm(Support ~ MNC1Final0 + MNC0Final1 + MNC1Final1 + CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                NetPPE + MarketCap + HighManufEnergy + electricpowergeneration, data = dat, family = 'binomial')

## Calculate pred probabilities for each case

holder <- as.data.frame(matrix(NA, 4, 3))
names(holder) <- c('Var', 'PredProb', 'SE')
holder[1:4, 1] <- c('MNC1Final0', 'MNC0Final1', 'MNC1Final1', 'MNC0Final0')

# case 1
newdat <- dat
newdat$MNC1Final0 <- 1
newdat$MNC0Final1 <- 0
newdat$MNC1Final1 <- 0

preds1 <- predict(fit1, newdat, type="response", se.fit=TRUE)
holder[1, 2] <- mean(preds1$fit)
holder[1, 3] <- sqrt(sum(preds1$fit^2))/length(preds1$se.fit)

# case 2
newdat$MNC1Final0 <- 0
newdat$MNC0Final1 <- 1
newdat$MNC1Final1 <- 0

preds2 <- predict(fit1, newdat, type="response", se.fit=TRUE)
holder[2, 2] <- mean(preds2$fit)
holder[2, 3] <- sqrt(sum(preds2$fit^2))/length(preds2$se.fit)

# case 3
newdat$MNC1Final0 <- 0
newdat$MNC0Final1 <- 0
newdat$MNC1Final1 <- 1

preds3 <- predict(fit1, newdat, type="response", se.fit=TRUE)
holder[3, 2] <- mean(preds3$fit)
holder[3, 3] <- sqrt(sum(preds3$fit^2))/length(preds3$se.fit)

# case 4
newdat$MNC1Final0 <- 0
newdat$MNC0Final1 <- 0
newdat$MNC1Final1 <- 0

preds4 <- predict(fit1, newdat, type="response", se.fit=TRUE)
holder[4, 2] <- mean(preds4$fit)
holder[4, 3] <- sqrt(sum(preds4$fit^2))/length(preds4$se.fit)

holder$Int <- holder$SE*1.96

holder$MNC <- c('MNC', 'Non-MNC', 'MNC', 'Non-MNC')
holder$Final <- c('Other', 'Final', 'Final', 'Other')

holder$Var <- factor(holder$Var, levels = c('MNC0Final0', 'MNC0Final1', 'MNC1Final0', 'MNC1Final1'))

# set colors
myblue <- "#4575B4"
myred <- "#D73027"

# plot results
ggplot(holder, aes(Var, PredProb, col = Final)) + geom_point() + geom_errorbar(aes(Var, ymin = PredProb - Int, ymax = PredProb + Int), size = .75, width = .25) +
  labs(x = 'Non-MNC                                                   MNC   ', y = 'Average Predicted Probability') + 
  theme(axis.text.x=element_blank(),axis.ticks=element_blank(), legend.title= element_blank(), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
        panel.background = element_rect(fill = "white",colour = "black", size = 0.5, linetype = "solid"))  +
  scale_colour_manual(values = c(myblue, myred))

# Hypothesis test

t.test(preds1$fit,preds3$fit)

#######################################################################################
## APPENDIX B1
#######################################################################################
#######################################################################################
## FIGURE 10 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

statedata <- unique(fulldata[, c("State", "StateName", "CoalReserves")])
statedata$CoalReserves <- ifelse(is.na(statedata$CoalReserves), 0, statedata$CoalReserves)

names(statedata) <- c("state", "region", "Coal Reserves")
statedata$region <- as.character(statedata$region)
statedata$region <- tolower(statedata$region)
breaks <- quantile(statedata$`Coal Reserves`, probs = c(.2, .4, .6, .8, 1))


states <- map_data("state")
tfmerged <- merge(statedata, states, sort = FALSE, by = "region")
tfmerged <- tfmerged[order(tfmerged$order), ]
plot <- qplot(long, lat, data = tfmerged, group = group, fill = -log(`Coal Reserves` + 1), 
              geom="polygon", ylab = "", xlab = "")
plot + theme_bw() + theme(
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks = element_blank()) +
  labs(x="Coal Reserves by State", y="") + 
  theme(legend.position='none', axis.title.x = element_text(margin = margin(t = 16, r = 0, b = 0, l = 0)), text = element_text(size=20))




#######################################################################################
## FIGURE 11 
#######################################################################################


states <- map_data("state")
unique(states$region)
dat <- as.data.frame(matrix(c('connecticut', NA, NA,
                              'maine', NA, NA,
                              'massachusetts', NA, NA,
                              'new hampshire', NA, NA, 
                              'rhode island', NA, NA, 
                              'vermont', NA, NA, 
                              'new jersey', NA, NA, 
                              'new york', 737, 98.10,
                              'pennsylvania', 2230, 68.14,
                              'illinois', 2895, 47.47,
                              'indiana', 4697, 66.48,
                              'michigan', 891, 103.17,
                              'ohio', 1460, 89.59,
                              'wisconsin', 1519, 79.21,
                              'iowa', 2682, 52.04,
                              'kansas', NA, NA,
                              'minnesota', 1167, 53.24,
                              'missouri', 787, 59.55,
                              'nebraska', 392, 42.13,
                              'north dakota', NA, NA,
                              'south dakota', NA, NA,
                              'delaware', NA, NA,
                              'district of columbia', NA, NA,
                              'florida', 933, 97.69,
                              'georgia', 1045, 104.97,
                              'maryland', 909, 58.43,
                              'north carolina', 869, 106.64,
                              'south carolina', 896, 95.81,
                              'virginia', 1641, 99.77,
                              'west virginia', 764, 106.03,
                              'alabama', 1341, 91.43,
                              'kentucky', 1026, 88.90,
                              'mississippi', NA, NA,
                              'tennessee', 2525, 92.24,
                              'arkansas', 298, 89.39,
                              'louisiana', NA, NA,
                              'oklahoma', NA, NA,
                              'texas', 833, 79.84,
                              'arizona', 431, 52.59,
                              'colorado', NA, NA,
                              'idaho', 414, 50.69,
                              'montana', NA, NA,
                              'nevada', NA, NA,
                              'new mexico', NA, NA,
                              'utah', 718, 54.58,
                              'wyoming', 1553, 31.21,
                              'alaska', 4, NA,
                              'california', 1330, 69.39,
                              'hawaii', NA, NA,
                              'oregon', NA, NA,
                              'washington', NA, NA), ncol = 3, byrow = T))

names(dat) <- c("region", "IndustrialCoalConsumption", "IndustrialCoalPrice")
dat$region <- as.character(dat$region)
dat <- dat[order(dat$region), ]

gdp <- c(138090,
         37011,
         215335,
         84227,
         1623600,
         218481,
         203970,
         54889,
         64827,
         639573,
         338522,
         49808,
         45733,
         568001,
         232204,
         119818,
         103760,
         130354,
         182018,
         42762,
         233004,
         326902,
         323117,
         230313,
         76758,
         206526,
         29069,
         74767,
         111244,
         52736,
         418140,
         61562,
         973550,
         346777,
         27123,
         407741,
         115657,
         145716,
         490192,
         41137,
         130248,
         33525,
         213822,
         1008630,
         95406,
         21038,
         333522,
         280758,
         50037,
         212779,
         31747)

dat <- cbind(dat, gdp)
dat$IndustrialCoalConsumption <- as.numeric(as.character(dat$IndustrialCoalConsumption))
dat$ScaledCoalConsumption <- dat$IndustrialCoalConsumption / dat$gdp

dat <- merge(dat, states, sort = FALSE, by = 'region')
dat <- dat[order(dat$order), ]
dat$IndustrialCoalConsumption <- as.numeric(as.character(dat$IndustrialCoalConsumption))
dat$IndustrialCoalPrice <- as.numeric(as.character(dat$IndustrialCoalPrice))
dat$ScaledCoalConsumption <- as.numeric(as.character(dat$ScaledCoalConsumption))

plot <- qplot(long, lat, data = dat, group = group, fill = -log(IndustrialCoalConsumption), 
              geom="polygon", ylab = "", xlab = "")
plot1 <- plot + theme_bw() + theme(
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks = element_blank()) +
  labs(x="Industrial Coal Consumption (Scaled by GDP)", y="") + 
  theme(legend.position='none', axis.title.x = element_text(margin = margin(t = 16, r = 0, b = 0, l = 0)), text = element_text(size=20))


plot <- qplot(long, lat, data = dat, group = group, fill = -log(IndustrialCoalPrice), 
              geom="polygon", ylab = "", xlab = "")
plot2 <- plot + theme_bw() + theme(
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks = element_blank()) +
  labs(x="Average Coal Price (Industrial End Use)", y="") + 
  theme(legend.position='none', axis.title.x = element_text(margin = margin(t = 16, r = 0, b = 0, l = 0)), text = element_text(size=20))

grid.arrange(plot1, plot2, ncol = 2)

#######################################################################################
## FIGURE 12 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Plot lobbying positions
fulldata_sub <- na.omit(fulldata[, c("HR2454_Pos", "Industry")])
fulldata_sub$HR2454_Pos <- ifelse(fulldata_sub$HR2454_Pos == 1, "For", ifelse(fulldata_sub$HR2454_Pos == 0, "Against", fulldata_sub$HR2454_Pos))
fulldata_sub$HR2454_Pos <- factor(fulldata_sub$HR2454_Pos, levels = c("For", "Against", "99"))
pro <- fulldata_sub[fulldata_sub$HR2454_Pos == "For" | fulldata_sub$HR2454_Pos == "Against", ]
pro$dummy <- 1
countsfor <- tapply(pro$dummy, pro$Industry, sum)
countsfor <- sort(countsfor, decreasing = T)
pro$Industry <- factor(pro$Industry, levels = names(countsfor))

ggplot(pro, aes(Industry)) +
  geom_bar(aes(fill = HR2454_Pos), ) + 
  labs(y="Firms Lobbying on H.R.2454") + scale_fill_manual("legend", values=c("#4575B4", "#D73027")) + 
  coord_flip() + theme(legend.title=element_blank())



#######################################################################################
## TABLE 1 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Model 1
fit1 <- glm.cluster(data = fulldata, Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results1 <- summary(fit1)[, c(1,2,4)]

# Model 2
fit2 <- glm.cluster(data = fulldata, Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results2 <- summary(fit2)[, c(1,2,4)]

# Assessing substantive magnitude for model 1
fakedata <- fulldata
fakedata$CompShareZipCarbon <- mean(fakedata$CompShareZipCarbon, na.rm = T)
fakedata2 <- fulldata
fakedata2$CompShareZipCarbon <- mean(fulldata$CompShareZipCarbon, na.rm = T) + sd(fulldata$CompShareZipCarbon, na.rm = T)
fakeprob <- predict.glm(glm(Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                              NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration,
                            data = fulldata, family = "binomial"), fakedata, type = 'response')
fakeprob2 <- predict.glm(glm(Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                               NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration,
                             data = fulldata, family = "binomial"), fakedata2, type = 'response')
change <- mean(fakeprob2 - fakeprob, na.rm = T)
print(change / mean(fulldata$Support))

# Assessing substantive magnitude for model 2
fakedata <- fulldata
fakedata$CompShareZipCoal <- mean(fakedata$CompShareZipCoal, na.rm = T)
fakedata2 <- fulldata
fakedata2$CompShareZipCoal <- mean(fulldata$CompShareZipCoal, na.rm = T) + sd(fulldata$CompShareZipCoal, na.rm = T)
fakeprob <- predict.glm(glm(Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                              NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, data = fulldata, family = "binomial"), fakedata, type = 'response')
fakeprob2 <- predict.glm(glm(Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                               NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, data = fulldata, family = "binomial"), fakedata2, type = 'response')
change <- mean(fakeprob2 - fakeprob, na.rm = T)
print(change / mean(fulldata$Support))

# Print table
table_me(list(results1, results2), list(nrow(fit1[[1]]$model), nrow(fit2[[1]]$model)), "Baseline Results", 2)

#######################################################################################
## TABLE 2 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Model 1
fit1 <- zelig(Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, 
              data = fulldata, model = "relogit")

results1 <- cbind(unlist(fit1$get_coef()), unlist(fit1$get_se()), unlist(fit1$get_pvalue()))

# Model 2
fit2 <- zelig(Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, 
              data = fulldata, model = "relogit")

results2 <- cbind(unlist(fit2$get_coef()), unlist(fit2$get_se()), unlist(fit2$get_pvalue()))

# Excluding outliers

# Model 3
dummy.fit1 <- glm(Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + NetPPE + MarketCap + HighManufEnergy + FSubs + 
              electricpowergeneration, data = fulldata, family = "binomial"(link = "logit"))

resid <- residuals(dummy.fit1)
threetimes.a <- sd(resid) * 3

fit3 <- glm.cluster(data = fulldata[resid < threetimes.a, ], formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results3 <- summary(fit3)[, c(1,2,4)]

# Model 4
dummy.fit2 <- glm(Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + NetPPE + MarketCap + HighManufEnergy + FSubs + 
              electricpowergeneration, data = fulldata, family = "binomial"(link = "logit"))

resid <- residuals(dummy.fit2)
threetimes.b <- sd(resid) * 3

fit4 <- glm.cluster(data = fulldata[resid < threetimes.b, ], formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results4 <- summary(fit4)[, c(1,2,4)]

table_me(list(results1, results2, results3, results4), 
         list(length(fit1$get_residuals()[[1]]), length(fit2$get_residuals()[[1]]), 
              nrow(fit3[[1]]$model), nrow(fit4[[1]]$model)), 'Additional Results', 2)

#######################################################################################
## TABLE 3 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# model 1
fit1 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + allocation, cluster = "CAICS12COD", family = "binomial")

results1 <- summary(fit1)[, c(1,2,4)]

# model 2
fit2 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + allocation, cluster = "CAICS12COD", family = "binomial")

results2 <- summary(fit2)[, c(1,2,4)]
        
# model 3
fit3 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + FSubs + electricpowergeneration + teii, cluster = "CAICS12COD", family = "binomial")

results3 <- summary(fit3)[, c(1,2,4)]

# model 4
fit4 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + FSubs + electricpowergeneration + teii, cluster = "CAICS12COD", family = "binomial")

results4 <- summary(fit4)[, c(1,2,4)]

# model 5
fit5 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + electricalmanufacturing, cluster = "CAICS12COD", family = "binomial")

results5 <- summary(fit5)[, c(1,2,4)]

# model 6
fit6 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + electricalmanufacturing, cluster = "CAICS12COD", family = "binomial")

results6 <- summary(fit6)[, c(1,2,4)]

table_me(list(results1, results2, results3, results4, results5, results6), 
         list(nrow(fit1[[1]]$model), nrow(fit2[[1]]$model), 
              nrow(fit3[[1]]$model), nrow(fit4[[1]]$model),
              nrow(fit5[[1]]$model), nrow(fit6[[1]]$model)), 'Lobbying for Private Goods', 2)


#######################################################################################
## TABLE 4 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Model 1 
fit1 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + RES, cluster = "CAICS12COD", family = "binomial")

results1 <- summary(fit1)[, c(1,2,4)]

# Model 2
fit2 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + RES, cluster = "CAICS12COD", family = "binomial")

results2 <- summary(fit2)[, c(1,2,4)]

# Model 3 
fit3 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + RES.Years, cluster = "CAICS12COD", family = "binomial")

results3 <- summary(fit3)[, c(1,2,4)]

# Model 4
fit4 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + RES.Years, cluster = "CAICS12COD", family = "binomial")

results4 <- summary(fit4)[, c(1,2,4)]

# Model 5
fit5 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + RES*NetPPE, cluster = "CAICS12COD", family = "binomial")

results5 <- summary(fit5)[, c(1,2,4)]

# Model 6
fit6 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + RES*NetPPE, cluster = "CAICS12COD", family = "binomial")

results6 <- summary(fit6)[, c(1,2,4)]

table_me(list(results1, results2, results3, results4, results5, results6), 
         list(nrow(fit1[[1]]$model), nrow(fit2[[1]]$model), 
              nrow(fit3[[1]]$model), nrow(fit4[[1]]$model),
              nrow(fit5[[1]]$model), nrow(fit6[[1]]$model)), 'Dynamics of Climate Change Policy', 2)

#######################################################################################
## TABLE 5 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Model 1
fit1 <- glm.cluster(data = fulldata, formula = Support ~ TotalAssets*CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results1 <- summary(fit1)[, c(1,2,4)]

# Model 2
fit2 <- glm.cluster(data = fulldata, formula = Support ~ TotalAssets*CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results2 <- summary(fit2)[, c(1,2,4)]

# Model 3
fit3 <- glm.cluster(data = fulldata, formula = Support ~ Productivity*CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results3 <- summary(fit3)[, c(1,2,4)]

# Model 4
fit4 <- glm.cluster(data = fulldata, formula = Support ~ Productivity*CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results4 <- summary(fit4)[, c(1,2,4)]

table_me(list(results1, results2, results3, results4), 
         list(nrow(fit1[[1]]$model), nrow(fit2[[1]]$model), 
              nrow(fit3[[1]]$model), nrow(fit4[[1]]$model)), 'Interaction with Firm Size and Productivity', 2)


#######################################################################################
## TABLE 6 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Model 1
fit1 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal*LogPenetration + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results1 <- summary(fit1)[, c(1,2,4)]

# Model 2
fit2 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal*teii + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results2 <- summary(fit2)[, c(1,2,4)]

# Model 3
fit3 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal*LogExports + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results3 <- summary(fit3)[, c(1,2,4)]

# Model 4
fit4 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon*LogExports + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results4 <- summary(fit4)[, c(1,2,4)]

# Model 5
fit5 <- brglm(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs*FinalGood + electricpowergeneration, cluster = "CAICS12COD", family = "binomial")

results5 <- summary(fit5)$coefficients[, c(1,2,4)]

table_me(list(results1, results2, results3, results4, results5), 
         list(nrow(fit1[[1]]$model), nrow(fit2[[1]]$model), 
              nrow(fit3[[1]]$model), nrow(fit4[[1]]$model),
              nrow(fit5$model)), 'Interaction with Trade Exposure and Production Stage Characteristics', 2)


#######################################################################################
## APPENDIX B3
#######################################################################################

#######################################################################################
## Figure 14 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

dat <- fulldata %>%
  group_by(CAICS12COD) %>%
  summarize(sd.compshare = sd(CompShareZipCarbon, na.rm = T), mean.compshare = mean(CompShareZipCarbon, na.rm = T),
            mean.assets = mean(TotalAssets, na.rm = T), n.firms = n(), support = max(Support, na.rm = T),
            mean.ppe = mean(NetPPE, na.rm = T), mean.zipenergy = mean(HighZipEnergy, na.rm = T), sd.zipenergy = sd(HighZipEnergy, na.rm = T))

dat$support <- ifelse(dat$support == 1, "Support", "No Support")
dat$support <- factor(dat$support, levels = c("Support", "No Support"))

#### industry level
# assets
plot1 <- ggplot(dat, aes(log(mean.assets), mean.zipenergy, col = support)) +
  geom_point() + theme(legend.title=element_blank()) +
  labs(x="(Log) Mean Total Assets", y = "Prop. High Cost") +
  theme(text = element_text(size=25)) 

# ppe
plot2 <- ggplot(dat, aes(mean.ppe, mean.zipenergy, col = support)) +
  geom_point() + theme(legend.title=element_blank()) +
  labs(x="Mean PP&E", y = "Prop. High Cost") +
  theme(text = element_text(size=25)) 

# n firms
plot3 <- ggplot(dat, aes(n.firms, mean.zipenergy, col = support)) +
  geom_point() + theme(legend.title=element_blank()) +
  labs(x="No. of Firms", y = "Prop. High Cost") +
  scale_x_continuous(limits = c(0, 400)) +
  theme(text = element_text(size=25)) 

grid.arrange(plot1, plot2, plot3, ncol = 1)

#######################################################################################
## Figure 15 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Top panel
plot1 <- ggplot(fulldata, aes(log(TotalAssets), HighZipEnergy)) +
  geom_point() +
  labs(x="(Log) Total Assets", y = "High Costs") +
  theme(text = element_text(size=25)) 

# Bottom panel
plot2 <- ggplot(fulldata, aes(NetPPE, HighZipEnergy)) +
  geom_point() +
  labs(x="Net PP&E", y = "High Costs") +
  scale_x_continuous(limits = c(0, 40000000)) +
  theme(text = element_text(size=25)) 

grid.arrange(plot1, plot2, ncol = 1)

#######################################################################################
## Figure 16
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Top panel
fit1 <- logitmfx(Support ~ CompShareHS + HighEnergy + MarketShare + Productivity + NetPPE + 
                   MarketCap + HighManufEnergy + FSubs + electricpowergeneration, 
                 data = fulldata, atmean = T, clustervar1 = "CAICS12COD")
mfx <- as.data.frame(fit1$mfxest)

z_score <- qnorm(0.975)
mfx$UpperCI <- mfx[,1] + mfx[,2] * z_score
mfx$LowerCI <- mfx[,1] - mfx[,2] * z_score
mfx$Variable <- row.names(mfx)
names(mfx) <- c("MarginalEffect", "StError", "z", "P", "UpperCI", "LowerCI", "Variable")

# Re-name variables
mfx[1,7] <- "Competitor Costs (State-Level Carbon Reliance)"
mfx[2,7] <- "Own Costs (State-Level Carbon Reliance)"
mfx[5,7] <- "PPE"
mfx[7,7] <- "Energy Intensive Manufacturing"
mfx[8,7] <- "Multinational"
mfx[9,7] <- "Electricity Provider"

# Re-order factors
mfx[,7] <- factor(mfx[,7], levels = c(
  "PPE",
  "MarketCap",
  "Productivity",
  "MarketShare",
  "Electricity Provider",
  "Energy Intensive Manufacturing",
  "Multinational",
  "Own Costs (State-Level Carbon Reliance)",
  "Competitor Costs (State-Level Carbon Reliance)"
))

# Set colors
myblue <- "#4575B4"


# coefplot
stateplot1 <- ggplot(mfx[c(1,2,8, 3),], aes(Variable, MarginalEffect)) +
  scale_x_discrete('') +
  # scale_y_continuous('', limits = c(-.05, .08)) +
  theme_bw() + 
  geom_errorbar(aes(x = Variable, ymin = LowerCI, ymax = UpperCI), size = 2, width = .5, color = myblue) + 
  geom_point(aes(x = Variable, y = MarginalEffect), size = 5, color = myblue) +
  geom_hline(yintercept = 0) + 
  theme(text = element_text(size=20)) +
  coord_flip() 

# Bottom panel

fit1 <- logitmfx(Support ~ CompShareStateCoal + CoalReserves + MarketShare + Productivity + NetPPE + 
                   MarketCap + HighManufEnergy + FSubs + electricpowergeneration, 
                 data = fulldata, atmean = T, clustervar1 = "CAICS12COD")
mfx <- as.data.frame(fit1$mfxest)

z_score <- qnorm(0.975)
mfx$UpperCI <- mfx[,1] + mfx[,2] * z_score
mfx$LowerCI <- mfx[,1] - mfx[,2] * z_score
mfx$Variable <- row.names(mfx)
names(mfx) <- c("MarginalEffect", "StError", "z", "P", "UpperCI", "LowerCI", "Variable")

# Re-name variables
mfx[1,7] <- "Competitor Costs (State-Level Coal Reserves)"
mfx[2,7] <- "Own Costs (State-Level Coal Reserves)"
mfx[5,7] <- "PPE"
mfx[7,7] <- "Energy Intensive Manufacturing"
mfx[8,7] <- "Multinational"
mfx[9,7] <- "ElectricityProvider"

# Re-order factors
mfx[,7] <- factor(mfx[,7], levels = c(
  "PPE",
  "MarketCap",
  "Productivity",
  "MarketShare",
  "ElectricityProvider",
  "Energy Intensive Manufacturing",
  "Multinational",
  "Own Costs (State-Level Coal Reserves)",
  "Competitor Costs (State-Level Coal Reserves)"
))


# Set colors
myred <- "#D73027"


# coefplot
stateplot2 <- ggplot(mfx[c(1,2,8,3),], aes(Variable, MarginalEffect)) +
  scale_x_discrete('') +
  theme_bw() + 
  geom_errorbar(aes(x = Variable, ymin = LowerCI, ymax = UpperCI), size = 2, width = .5, color = myred) + 
  geom_point(aes(x = Variable, y = MarginalEffect), size = 5, color = myred) +
  geom_hline(yintercept = 0) + 
  theme(text = element_text(size=20)) +
  coord_flip() 

grid.arrange(stateplot1, stateplot2, ncol = 1)



#######################################################################################
## Figure 17 
#######################################################################################

set.seed(123)

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

##### (Fake) competitor vars using carbon intensity (randomize NetGenerationFromCoal)
vars <- c('Support', 'CompShareZipCarbon', 'HighZipEnergy', 'MarketShare', 'Productivity',
          'NetPPE', 'MarketCap', 'HighManufEnergy', 'FSubs', 'electricpowergeneration')

holder <- NA

for(j in 1:1000){
  dat <- fulldata[complete.cases(fulldata[, vars]), ]
  
  N <- round(nrow(dat) * .15, digits = 0)
  firms <- sample(1:nrow(dat), N, replace = F)
  for(i in 1:length(firms)){
    dat[firms[i], 'NetGenerationFromCoal'] <- sample(dat$NetGenerationFromCoal, 1)
    
  }
  
  dat$Fake_CompShareZipCarbon <- NA
  dat$Fake_HighZipEnergy <- ifelse(dat$NetGenerationFromCoal > median(dat$NetGenerationFromCoal, na.rm = T), 1, 0)
  
  for(i in 1:nrow(dat)){
    ind <- dat[i,"CAICS12COD"]
    temp <- dat[fulldata$NAME != dat$NAME[i] & dat$CAICS12COD == ind, ]
    high <- temp[temp$Fake_HighZipEnergy == 1, ]
    
    # prop of competitors in high energy states
    dat$Fake_CompShareZipCarbon[i] <- length(na.omit(high$NAME)) / length(na.omit(temp$NAME))
    
  }
  
  fake_fit <- summary(glm.cluster(data = dat, formula = Support ~ Fake_CompShareZipCarbon + Fake_HighZipEnergy + MarketShare + Productivity + 
                                    NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial"))
  
  holder <- c(holder, fake_fit[2,1])
}

dist1 <- na.omit(holder)
quant1 <- quantile(dist1, probs = .99)

real_fit <- summary(glm.cluster(data = dat, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                                  NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial"))

real1 <- real_fit[2,1]

#### (Fake) Competitor vars using coal production (randomize ZipCoalField)

vars <- c('Support', 'CompShareZipCoal', 'ZipCoalField', 'MarketShare', 'Productivity',
          'NetPPE', 'MarketCap', 'HighManufEnergy', 'FSubs', 'electricpowergeneration')

holder <- NA

for(j in 1:1000){
  dat <- fulldata[complete.cases(fulldata[, vars]), ]
  N <- round(nrow(dat) * .15, digits = 0)
  firms <- sample(1:nrow(dat), N, replace = F)
  dat$Fake_ZipCoalField <- dat$ZipCoalField
  
  for(i in 1:length(firms)){
    dat[firms[i], 'Fake_ZipCoalField'] <- sample(dat$ZipCoalField, 1)
    
  }
  
  dat$Fake_CompShareZipCoal <- NA
  
  # Loop across firms
  
  for(i in 1:nrow(dat)){
    ind <- dat[i,"CAICS12COD"]
    temp <- dat[dat$NAME != dat$NAME[i] & dat$CAICS12COD == ind, ]
    
    high <- temp[temp$Fake_ZipCoalField == 1, ]
    dat$Fake_CompShareZipCoal[i] <- length(na.omit(high$NAME)) / length(na.omit(temp$NAME))
    
  }
  
  
  fake_fit <- summary(glm.cluster(data = dat, formula = Support ~ Fake_CompShareZipCoal + Fake_ZipCoalField + MarketShare + Productivity + 
                                    NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial"))
  
  holder <- c(holder, fake_fit[2,1])
}

dist2 <- na.omit(holder)
quant2 <- quantile(dist2, probs = .99)

real_fit <- summary(glm.cluster(data = dat, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                                  NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration, cluster = "CAICS12COD", family = "binomial"))

real2 <- real_fit[2,1]

perm <- as.data.frame(cbind(dist1, dist2))

# Left panel
plot1 <- ggplot(perm, aes(dist1)) + geom_density() + 
  geom_vline(xintercept = quant1, col = 'red', lty = 2) +
  geom_vline(xintercept = real1, col = 'red') + 
  labs(x = 'Coefficient Estimate Competitor Costs (Coal Reliance)', y = 'Density') + 
  theme(text = element_text(size=15)) 



plot2 <- ggplot(perm, aes(dist2)) + geom_density() + 
  geom_vline(xintercept = quant2, col = 'red', lty = 2) +
  geom_vline(xintercept = real2, col = 'red') + 
  labs(x = 'Coefficient Estimate Competitor Costs (Coal Reserves)', y = 'Density') + 
  theme(text = element_text(size=15)) 

grid.arrange(plot1, plot2, ncol = 2)


#######################################################################################
## Table 7 
#######################################################################################

# Re-load main data file
fulldata <- read.csv('Enemies_Replication.csv')

# Model 1
fit1 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + USNWR_BusinessIndex, cluster = "CAICS12COD", family = "binomial")

results1 <- summary(fit1)[, c(1,2,4)]

# Model 2
fit2 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + USNWR_BusinessIndex, cluster = "CAICS12COD", family = "binomial")

results2 <- summary(fit2)[, c(1,2,4)]

# Model 3
fit3 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + USNWR_TaxBurden, cluster = "CAICS12COD", family = "binomial")

results3 <- summary(fit3)[, c(1,2,4)]

# Model 4
fit4 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + USNWR_TaxBurden, cluster = "CAICS12COD", family = "binomial")

results4 <- summary(fit4)[, c(1,2,4)]

# Model 5
fit5 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCarbon + HighZipEnergy + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + Delaware, cluster = "CAICS12COD", family = "binomial")

results5 <- summary(fit5)[, c(1,2,4)]

# Model 6
fit6 <- glm.cluster(data = fulldata, formula = Support ~ CompShareZipCoal + ZipCoalField + MarketShare + Productivity + 
                      NetPPE + MarketCap + HighManufEnergy + FSubs + electricpowergeneration + Delaware, cluster = "CAICS12COD", family = "binomial")

results6 <- summary(fit6)[, c(1,2,4)]

table_me(list(results1, results2, results3, results4, results5, results6), 
         list(nrow(fit1[[1]]$model), nrow(fit2[[1]]$model), 
              nrow(fit3[[1]]$model), nrow(fit4[[1]]$model),
              nrow(fit5[[1]]$model), nrow(fit6[[1]]$model)), 'Headquarter Incentives', 2)


